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ABSTRACT 

This report outlines the development of a general purpose aerodynamic solver for 
compressible turbulent flows. Turbulent closure is achieved using either two-equation 
or Reynolds stress transport equations. The applicable equation set consists of Favre- 
averaged conservation equations for the mass, momentum and total energy, and trans- 
port equations for the turbulent stresses and turbulent dissipation rate. In order to 
develop a scheme with good shock-capturing capabilities, good accuracy and general 
geometric capabilities, a multi-block, cell-centered finite-volume approach is used. A 
Roe flux-difference splitting technique, coupled with a MUSCL scheme, is applied to 
the system of conservation and transport equations. Viscous fluxes are discretized 
using a finite- volume representation of a central-difference operator and the source 
terms are treated as an integral over the control-volume. The methodology is vali- 
dated by testing the algorithm on both two-dimensional and three-dimensional flows. 
Both the two-equation and Reynolds stress models are used on a two-dimensional 
10° compression ramp at Mach 3, and the two-equation model is used on the three- 
dimensional flow over a cone at angle of attack at Mach 3.5. With the development of 
this algorithm, it is now possible to compute complex, compressible high-speed flow 
fields using both two-equation and Reynolds stress turbulent closure models, with the 
capability of eventually evaluating their predictive performance. 
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1 Introduction 

There has been an increasing desire to predict the turbulent flow behavior over 
complicated aerodynamic geometries in both low- and high-speed flows. This goal 
will certainly lie outside the scope of direct numerical simulations (DNS) and even 
large eddy simulations (LES) in the near future [1], even with anticipated increases 
in computational speed. Thus practical computations will necessarily rely on the 
utilization of phenomenological models for the turbulent transport equations. Until 
recently, work in these areas was further limited to the incompressible regime with 
relatively little focused research in the compressible area [2]. In the incompressible 
regime, there has been significant progress in the area of predictive solutions for tur- 
bulent flows in a wide variety of cases to date. In the more complex flow situations, 
the simpler turbulence models can be readily implemented and these have been used 
most extensively. Their success has been measured by their performance in generat- 
ing mean flows and wall pressures, which are not strongly dependent on the correct 
modeling of turbulent transport properties. However, the complex flows where tur- 
bulent transport is important require more sophisticated turbulence models in order 
to capture the underlying physics. 

It is only with the initiation of such national priority efforts as the National 
Aero-Space Plane (NASP) and the High-Speed Civil Transport (HSCT) that there 
has been a strong impetus into the compressible regime for these turbulent stress 
closure models. Unfortunately, the practical aerodynamic flows that we are ultimately 
interested in are three dimensional with strong vortical regions and generally far from 
equilibrium. As has already been noted, such fully three-dimensional flows cannot be 
adequately predicted using simple eddy viscosity type models because of the strong 
anisotropies in the flow. It is necessary to increase the physics within the turbulent 
transport models, and this can be done by advancing to Reynolds stress transport 
models. While these models contain improved capabilities for the stress anisotropies 
and the associated strain histories, their computational stiffness has lessoned their 
appeal to the general user community. These stiffness problems can be traced to the 
near-wall corrections applied to the high Reynolds number versions of the transport 
models. For this reason, wall function corrections have proved useful since they can 
remove the stiffness problems; however, this is done at the expense of constraining 
the Reynolds stress transport models to attached flows. 

The task at hand is to develop a versatile means of implementing Reynolds stress 
transport models into generally applicable aerodynamic codes for application to high- 
speed compressible flows. This task, while straightforward in concept, is practically 
very cumbersome because of the significant number of terms in the turbulent trans- 
port equations and the nonlinear coupling between the mean flow equations and the 
transport equations. Add to this the unknown correlations resulting from the com- 
pressibilty of the flow and one can see that it is a challenging task to effectively 
develop a robust computational tool for studying complex compressible flows using 
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Reynolds stress turbulence models. 

Mean conservation and transport equations used in the solution of compressible 
turbulent flows are presented in Section 2. These include the conservation of mean 
mass, momentum and energy, as well as the transport equations for the various tur- 
bulent Reynolds stresses and turbulent dissipation rate. The models for the unknown 
turbulent correlations in the high Reynolds number form of these equations are also 
presented. The corrections to these models to account for the near-wall effects of 
solid boundaries is discussed separately in Section 3. In Section 4, the 12 partial 
differential equations (7 for the two-equation model) comprising the solution set are 
recast in vector form. The dependent variable vector includes the density, velocities, 
total energy, turbulent Reynolds stresses and dissipation rate. The discretization of 
the resulting inhomogeneous flux equation, expressed in generalized coordinates, is 
then discussed. 

In the context of the present report where the intent is to establish the differential 
and numerical framework for the solution o f th e equations, it is not advantageous to do 
a detailed predictive study of particular compressible flows. It is necessary, however, 
to compute some representative flow fields in order to validate the capabilities of 
the numerical code. Detailed comparisons will require the choice of a more accurate 
near-wall model for the Reynolds stress formulation. The two flows that will be 
considered here are the two-dimensional compression ramp and the three-dimensional 
cone. For the ramp problem, the flow is computed using both the two-equation and 
Reynolds stress models, and for the cone problem only the two-equation model is 
used. The results from the calculations are given in Section 5 where profiles at 
selected streamwise (ramp) or azimuthal (cone) stations are shown for both the mean 
and turbulent quantities. 
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2 Conservation and Transport Equations 

It has been known for some time that the optimal way, both mathematically and 
physically, of formulating both the mean conservation equations and the Reynolds 
stress transport equations is to employ Favre or mass- weighted averages [3, 4], This 
is due to the resulting similarity between the incompressible and compressible form 
of the equations. For completeness, the relevant mean conservation and Reynolds 
stress transport equations are given here. Any dependent flow variable / can be 
decomposed using the usual Reynolds decomposition given by 


/ = / + /', ( 2 . 1 ) 

where 

or a Favre-average given by 

/ = / + /", ( 2 . 2 ) 

where 


I 


pj_ 

p 


(2.2a) 


In the equation development that follows, the dependent variables are the density 
p, the pressure p, the temperature T, the velocity iq, the turbulent Reynolds stress 
Tij and the turbulent dissipation rate e ir The differential equations are written in 
Cartesian tensor notation for compactness and consistency. In the following deriva- 
tions, the sign for the unknown correlations means the exact functional form, 
and the sign means the modeled or approximate form. 


2.1 Conservation Equations 

The mean conservation equations for mass, momentum and total energy can then be 
written as 

MASS 

d t p + (puk),k = 0 (2.3) 


MOMENTUM 

9 t (pui) + (puiUj + pSij),j = aijj - ( pT tj ),j (2.4) 

where 

& ij ~ T T ^j f i) 
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2 

~ IUk,kfiij + 7 *(« ,j + Uj,i) (2.4a) 

o 

is the viscous stress tensor with p the mean molecular viscosity calculated from 
Sutherland’s law. 

TOTAL ENERGY 


d t (pE) + [uj 

(pE + p)],j 

= (WijUi + aiju'l - 9j),j - pE" u j)ij 

(2.5) 

where 






- Uiiii 

e = c„t+ 2 + ' 2 ' 

(2.5a) 



q } = -ItT, j 




~ -kT , j 

(2.56) 


p E"u” 

~ ~pu n i u ll i u , \ 

- C v pu"T " + pUiTij + 2 , 

(2.5c) 


qj is the mean heat flux vector, and k is the mean thermal conductivity. The energy 
equation is written in terms of the total energy because this formulation is necessary 
in order to effectively employ shock capturing techniques in the numerical solution 
algorithm [5]. 


2.2 Transport Equations 

Equation (2.4) contains the Favre- averaged Reynolds stress tensor, (r,-j = u"u"). 
In the utilization of a two-equation closure, the individual stress components are 
obtained by way of a Boussinesq approximation relating the turbulent stresses to 
mean velocity gradients, 

PTij O' ^pkStj - 2 JI t Sij - ( 2 . 6 ) 

with the eddy viscosity given by 

F, = C r p~, (2.6a) 

where S,j(= (uij + u Ji ,)/2) is the strain rate tensor, and C M a constant usually taken 
to be (0.09). This type of closure is clearly inaccurate in regions close to solid bound- 
aries, since it assumes an isotropic distribution of the stresses with respect to the 
kinetic energy. This is easily demonstrated by recalling that in the near-wall region 
asymptotic analyses [6] show that while the streamwise and spanwise normal stresses 
behave as the kinetic energy (i.e. as 0(y 2 ), where y is the distance normal to the 
wall), the transverse normal stress r yy goes as 0(y 4 ). A near-wall correction is clearly 
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needed and this will be the topic of discussion in the next section for both the two- 
equation and Reynolds stress models. Nevertheless, this two-equation approach is 
quite useful since it couples the turbulent kinetic energy equation with an appropri- 
ate scale equation so that there are only two transport equations which need to be 
solved. 

The Reynolds stress transport equation approach solves transport equations for 
the individual stress components as well as an appropriate scale equation. While 
the deficiencies with the use of the Boussinesq approximation are removed, they are 
replaced with the complicating factor of needing to solve six transport equations for 
the various stress components (the Reynolds stress tensor being symmetric for flows of 
interest in aerodynamic flight). The incorporation of the Reynolds stress formulation 
within a general Navier-Stokes solver is intended to ultimately serve as a means of 
comparison with the two-equation formulation. 

Since the turbulent kinetic energy equation is simply the contraction of the full 
Reynolds stress transport equation, it suffices to outline the modeling of the Reynolds 
stress transport equation and point out any differences that may result. 

REYNOLDS STRESS TRANSPORT 


9i{pT tJ ) + (ukpTij),k = p i} + nj- + ng + Mij + n v l3 + d\ 3 

- pdj (2.7) 

where the right hand side represents the turbulent stress transport produced by the 
tuibulent production, P \j , the deviatoric part of the pressure strain-rate correlation, 
the pressure dilatation, II the mass flux variation, Mij y the viscous diffusion, 
the turbulent diffusion, D 1 -, and the turbulent dissipation rate. These terms are 
given by 

Pij = ~ PTikU h k ~ pT jk Ui,k 

(2.7a) 

n J = P'Kj + P' u j,i ~ 

(2.76) 

ng — -^p' u ’ ktk &ij 

(2.7c) 

M{j = U‘ (&jk,k H“ O'ik^ ~~ P)j —Pn ) 

(2.7(7) 

D ", = + <*>!),* 

(2.7c) 

D h = -[pu>"K + (p'u'Jjk + p'u' 3 8 lk )], k 

(2.7/) 

e V = OikUj'k + a 'jk u ',,k- 

(2.7s) 


The general form of the turbulent kinetic energy equation is obtained from Eq. (2.7) 
using the definition k = r,-,/2. 

The remaining transport equation that is used in both two-equation and Reynolds 
stress closures adopted here is the turbulent dissipation rate equation. Using the 
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defining relationship shown in Eq. (2.7g), it is derivable from the fluctuating mo- 
mentum equation. As in the incompressible case (e.g. [7]), this equation is quite 
complex and there is no direct experimental evidence about any of the terms in the 
equation. The approach that has been taken to date in the development of a dissipa- 
tion rate model that can be used in compressible flows is to partition the dissipation 
rate into a solenoidal (incompressible) component and a compressible component 
[8, 9]. Compressibility effects are then represented as asymptotic corrections to the 
incompressible (solenoidal) dissipation rate. This allows for the direct utilization of 
the incompressible form of the dissipation rate transport equation. The form of the 
dissipation rate partitioning that is adopted in the present work is given by [9] 

e = e s + e c 

— £s(l + c*i A/ t 2 ) (2-8) 

where a i( = 0.6) is a numerical constant extracted from comparisons with direct 
simulation data [2], and M t {— ^/EJ/floo) is the turbulent Mach number. A general 
high Reynolds number form of the isotropic solenoidal dissipation rate equation has 
been derived [10]. 

DISSIPATION RATE EQUATION 

4 1 D j/ 

dt{pts) + {S^kp^s\k = ~p£s&i,i + ~ ^ c 9 (2.9) 

where P Cs is the production, T) €s is the destruction, D[ s is the turbulent diffusion, and 
D v f is the viscous diffusion. 


2.3 Higher-Order Correlation Models 

As can be seen from an examination of the compressible mean conservation and tur- 
bulent transport equations, there are several higher order correlations which need to 
be modeled. In the incompressible regime, several of these terms have well estab- 
lished models which have been tested and proven in a variety of flow situations. The 
modeling of the unknown correlations in the compressible case is not as mature so 
there are more uncertainties in the predictive capabilities of these equations. Since 
this is a very active research area at the present time, the focus of the present report 
is to develop the numerical techniques for a rather general form of the equations. 
As new models are developed, they then can be easily implemented into the existing 
numerical structure. 

In the equation for the total energy, Eq. (2.5), both the viscous diffusion and the 
energy flux need to be modeled. The viscous diffusion term that appears in the total 
energy equation is the contraction of the corresponding term in the Reynolds stress 
transport equation, Eq. (2.7e). In the solution of the Reynolds stress transport equa- 
tions in incompressible high Reynolds number flows, this term is frequently neglected 
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since it is usually dominated by turbulent diffusion effects. When it is taken into 
account, it is usually treated simply in its gradient diffusion form and is written as 


a ik U 'j + a jk U < - + Tjk.i + Tikj). (2.10) 

The contraction of Eq. (2.10) leads to the appropriate modef for the viscous diffusion 
term in the total energy equation. 

The energy flux term, Eq. (2.5c), requires models for both the heat flux and 
triple velocity correlation terms. At present there has not been any consistent effort 
to develop alternatives to the gradient diffusion hypothesis for the heat flux term. 
For simplicity, the heat flux term is then modeled as 


~pu"T" ~ —R t f,j (2.11) 

where R t is the thermal eddy diffusivity which is given here as (cf. Eq. (2.6a)), 


pP 


= C U 


Pr t e s 


th 

Pr t ' 


( 2 . 12 ) 


As is seen by the form of Eq. (2.12), the thermal eddy diffusivity used here is simply 
the usual eddy viscosity divided by the turbulent Prandtl number Pr<(= 0.9). 

The remaining contribution in the energy flux term is the triple correlation or the 
turbulent diffusion term. This term is also present in the Reynolds stress equation in 
its uncontracted form (cf. Eq. (2.7f). This term has also been traditionally modeled 
with a gradient diffusion hypothesis by using a turbulent eddy viscosity. The form 
adopted here is given by 


puy'u’l 




pk‘ 


( 7~ij,k T Tjk,i T k,j) 


(2.13) 


where C s is a numerical constant which is assigned the value of 0.18. This is an 
isotropized version of the model used by Launder, Reece, and Rodi [11]. In their 
model, the coefficient C s was chosen to be 0.11, but in its present form the higher 
value used here is more appropriate. This model incorporates the functional form of 
the eddy viscosity relationship defined previously but with a different proportionality 
coefficient, /I t /<r*., with cr). = 0.75. Contraction of Eq. (2.13) will lead to the required 
functional form used in the energy flux model. 

The models required for closure of the total energy equation have now been identi- 
fied and it is necessary to examine the remaining unknown correlations in the Reynolds 
stress transport equations. The mass flux variation appearing in the Reynolds stress 
transport equation is represented by Eq. (2.7d). Consistent with the other correla- 
tions which are unique to the compressible formulation, the mass flux term is also 
modeled by invoking the gradient diffusion hypothesis; 


C^P 

P@ p 


P,i = 


Pt 

P 2 <?P 


P,i 


(2.14) 
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where cr p is a constant whose value is 0.5. 

The modeling of the pressure dilatation term has been the subject of analysis 
following the partitioning ideas invoked earlier for the dissipation rate term [12, 13]. 
The form proposed in [12] is used here and is given by 

P' u 'k,k — (a 2 Tij u\j + a 3 ) (2.15) 

where a 2 (= 0.6) and a 3 (= 0.2) are numerical constants calibrated by comparison 
with direct numerical simulations [12]. 

The only remaining model that needs to be determined is for the deviatoric part 
of the pressure-strain rate correlation, Eq. (2.7b). At present there has been no work 
directed toward the development of a compressible pressure-strain rate correlation 
model. The approach taken has been to use variable density extensions of the incom- 
pressible form of the model. The compressibility effects have then been isolated into 
the pressure dilatation term which has been solely derived for compressible flows. In 
light of this approach and the fact that there are a significant number of pressure- 
strain rate models presented in the literature, it is sufficient in the present context 
to present the functional form of a commonly used model and show its incorporation 
into the numerical algorithm. One of the most commonly used pressure-strain rate 
correlation models is the model of Launder, Reece and Rodi [11]. While it is only 
linear in the anisotropies of the Reynolds stress, it has been used extensively on a 
variety of flows. The high Reynolds number form of the equation is given by 


n$ — 4>iji + <t>ij 2 

4>ij\ = —Cipebij 

ta = - \ p ^ 

30C 2 — 2 . _ 2 ~ r x 

— />«(«.,; + uj,i - -Uk,k(>ij) 

55 o 


(2.16) 

(2.16a) 


(2.166) 


where C x and C 2 are numerical constants given by 3.0 and 0.4, respectively, and 


bij = ^( Tn - | kSij ) (2.17a) 

Dij = -pTikUk,j - pTjkUk,i (2.176) 

It should be emphasized that the pressure-strain rate correlation is the subject of 
extensive research and that other models have been proposed which should perform 
better than the above model (e.g. [14] ). However, comparing a variety of closure 
models is not the intent of the present report, but is a course of study being actively 
pursued and will be reported on later. 

The remaining equation that needs to be modeled is the solenoidal dissipation rate 
equation, Eq. (2.9). This equation has been examined recently [15, 10] for application 
to compressible flows. The models for the terms in Eq. (2.9) are given by 
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(2.18a) 


1 Dv 

p ~m = [1 " m h “ *)1<< 


Cel ^ P ^ 1 

(2.186) 

= C c2 /f 

(2.18c) 

A'. = C, 

(2.18<i) 

D u = (?«•,*)»* 

(2.18e) 


where C c i,C e 2 and C t are modeling coefficients which take the values 1.50,1.83 and 
0.15, respectively, 7 is the ratio of specific heats (= 1.4), and 0.7) is the exponent 
of a power law approximation for the mean molecular viscosity [16]. 

Equation (2.9) can then be rewritten as 


d t (pe s ) + {pu k e s ), k = -[- + m ( 7 - l)]p€ s u^ - C n ~pT tJ {u lo - ]rUk,kf>ij) 


- C ^P~t + C t ( ) + (pe Sik ), k 

,k 


) 


(2.19) 


This now completes the specification of the transport equations and closure models 
needed for the solution of turbulent flows. The models presented up to this point have 
been high Reynolds number models, that is, models applicable to flows away from 
walls. In the presence of solid boundaries, the desire is to integrate directly to the wall 
so that wall functions do not have to be implemented. While the use of wall functions 
reduces the stiffness of the equation set, it constrains the application of the models 
to attached flows. Integrating directly to the wall requires near-wall modifications 
to the models already presented. For clarity, these near-wall modifications will be 
discussed in the next section for both the two-equation and Reynolds stress models. 
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3 Near- Wall Corrections 


The equation sets presented in the previous section outlined the high Reynolds 
number form of the equations. In the presence of solid boundaries, these equations 
need to be modified to account for the low Reynolds number flow near the solid 
surface. This area has received considerable attention in the incompressible regime 
(e.g. [17, 18]) while in the compressible regime work is only beginning [10]. Due to 
the uncertainty in developing such near-wall models at this point, it suffices here to 
show the types of modifications that are incorporated into the high Reynolds number 
models of the previous section. 

As was alluded to in the previous section, the two-equation turbulence models 
suffer from inaccurate distribution of Reynolds stress values in the vicinity of walls. 
This inaccuracy is accounted for by introducing a wall damping function into the 
defining relationship for the eddy viscosity, Eq. (2.6a). The modified near-wall form 
is 


f*t — 

where the damping function / M is given by 


(3.1) 


u = 




(3.2) 


where y + is the distance normal to the wall in wall variables, and Re t is the turbulent 
Reynolds number defined as P/De s . Note that far from the wall the value of 
is delimited by unity to consistently merge with the high Reynolds number form of 
the models. This analytic representation for f fl is by no means unique; however, 
comparison with direct simulation results [19] in incompressible boundary-layer flow 
has shown this to be a reasonable representation of the correct near-wall behavior. 
Note that it also carries over to the subsequent implementation of the eddy viscosity 
in the other closure approximation of the previous section. This includes the triple 
correlation model used in Eq. (2.13) which now becomes 


pu'fu'lu'l 


2 , k\ 

q sJn \Tij t k T Tjk,i T 7"ifc,j)- 

o €<> 


(3.3) 


or, for the two-equation model 


- -—(*,* + Tik,i). (3.4) 

i <j k 

where <7* = 0.75. 

The remaining term that needs to be modified in the equations to account for the 
presence of solid boundaries is the destruction term, in the solenoidal disspation 
rate equation. This term is then rewritten as 
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(3.5) 


where [6] 



f2 = 


1 — exp 




(3.6) 


Note that in the absence of / 2 the destruction term increases without bound as the 
wall is approached. 

For the Reynolds stress closure, a near-wall correction needs to be used for the 
pressure-strain rate correlation. This has been the topic of research in wall-bounded 
incompressible flows for some time [18], but with mixed success. Nevertheless, since 
the intent here is not to validate or compare models, it suffices to implement a near- 
wall correction which typifies the form and structure of a near-wall closure that can 
be implemented in the present numerical formulation. The near-wall modification 
that is adopted here is a variable density extension of the Shima model [20]. Even 
though this model has deficiencies which will be pointed out shortly, it is readily 
amenable to implementation in complex flows because it lacks any dependency on 
the wall normals [21] which cause ambiguities in complex flow situations. The Shima 
model introduces a near-wall correction to the pressure-strain correlation given by 


n,'j — 0 ljX T <Pij 2 T &ijw (3-7) 

where 

tf,, ^ -Cfptb,, (3.8) 

with 

C* = C x + (1 - C,)f w 
f w = exp[-(0.015\/fcy/f') 4 ], 

fa 2 given by Eq. (2.16b), and C x = 3.0. The coordinate y is measured normal 
surface. The near-wall correction is given by 

1 2 

<f>ij2 — [c*( F ij ~^Pkk^ij) T + Ujj ~Uk,k^ij')]fw> (3-9) 

where a = 0.45 and j3 — 0.08. 

A near-wall correction to the solenoidal dissipation rate needs to be implemented. 
For consisitency, the model proposed by Shima [20] is again used in the present study. 
The solenoidal dissipation rate equation is easily modified to include the near-wall 
Shima correction by replacing C t \ with a C* t defined as 

C* = C e i(i + /„), (3-10) 


(3.8a) 
(3.86) 
to the 
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and adding 


where 


D'f = 

£3 


C o — e a - v 


) k 2k j Jw 

(3.11) 

d 2 k 

dxidxi 

(3.11a) 


to the right hand side of Eq. (2.19). With the introduction of the near- wall model, 
the closure coefficients C t \ , C& and C e are given by 1.35,1.80 and 0.15, respectively, 
and f w is given by Eq. (3.8b). Since the destruction term, in the Shima model 
does not include any near- wall damping effect (/2 = 1, cf. Eq. (3.5)), it is readily ap- 
parent that the destruction term in the solenoidal dissipation rate equation increases 
unbounded as the wall is approached. The effect of this deficiency becomes more 
pronounced as the grid resolution is increased and, therefore, needs to be corrected 
for quantitative comparisons. 

The true test of these near-wall corrections lie in their application to a variety of 
wall-bounded compressible flows. The intent in this work is to develop the numerical 
framework for application to such flows and future work will focus on the testing 
of different types of models which are better posed and more suited for comparison 
purposes. 


4 Numerical Implementation 


In the previous sections, the mean Navier-Stokes and turbulent transport equa- 
tions were developed using Cartesian tensor notation. This notation compacted the 
form of the equations so that the various models could be discussed with some gen- 
erality. Of course the component equations need to be discretized - for a total of 7 
equations for a two-equation model and 12 equations for the Reynolds stress model. 
It is most convenient to discuss the discretization of these equations by defining vector 
arrays for the dependent variables and analyzing a set of vector equations. In ad- 
dition, for notational simplicity, the overbars and tildes are dropped^ from the mean 
variables with the understanding that these are Favre- averaged variables. 

Expanding Eqs. (2.3), (2.4), (2.5) and the modeled forms of (2.7) and (2.9), the 
mean Navier-Stokes and turbulent transport equations can be written in vector form 
in an arbitrary coordinate system as: 


OQ d(F-F v ) d(G-G v ) d(H-H v ) 

dt dt + drj + dc 

where Q is the vector of dependent variables 


Q = 


Q 

j= r 


P 

pu 

pv 

pw 

pE 

pi~ xx 
P T yy 

P T ZZ 


pT X y 

P^XZ 


P T yz 

, pc 


(4.1) 


(4.1a) 


and F,G,H are the inviscid (convective) fluxes, F V ,G V ,H V are the viscous (diffu- 
sive) fluxes, and 5 represents the source terms due to production, destruction, and 
redistribution. 

An equation of state is required to complete the system of equations. The perfect 
gas equation of state is used in this study: 


p = (~/-l)[pE - ^ p(u 2 + v 2 + w 2 ) - pk ] (4.2) 

The presence of the turbulent kinetic energy term in the equation of state arises 
from the Favre-avcraging and creates a strong coupling between the mean equations 
and the normal Reynolds stress components for the Reynolds stress models or the 
turbulent kinetic energy equation for the two-equation models. 
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4.1 Discretization 

The semi-discrete, finite-volume form of Eq. (4.1) is written as: 

(dQldt) ijk + [(F-F v )V(lJ} 1+1/2, j,* [(E Fv)^ (,/ 

+[(G - G v )Vr]/ J} itj+1/ 2 <k - [(G- G v )Wr]/J] itj _ 1/2ik 
+[(H-H v )V(/J]i J.M-1/2 [(-^ 1/2 

-(7 S )i* = 0 

where the fluxes are defined at the interfaces of the computational cell bounding the 
cell-average value, Qijk- The definition of the numerical flux function approximating 
the interface flux determines the characteristics of the numerical scheme. The current 
emphasis is to develop a method capable of calculating high-speed flow about complex 
configurations. Problems of interest include separated flows, flows with large pressure 
gradients, and shock/boundary-layer interactions. This requires that we develop a 
scheme with good shock-capturing capabilities, good accuracy and general geometric 
capabilities. 

Upwind schemes have proven reliable for calculating flows with strong shocks. 
Additionally, Roe’s approximate Riemann solver [22] has been shown to be accurate in 
viscous flows [23]. Upwind schemes allow the dissipation required of shock-capturing 
schemes to be scaled according to each wave type recognized by the flux function, 
thereby minimizing the numerical dissipation. (See e.g, [24] for an excellent review of 
upwind differencing techniques.) The Roe approximate Riemann solver for the mean 
flow equations coupled with two-equation turbulence models was previously derived 
and reported in [5]. Section 4.2 derives the scheme for the mean flow equations 
coupled with a Reynolds stress turbulence model. 

Second-order spatial accuracy for the inviscid terms is attained by using the 
MUSCL scheme of van Leer [25]. The variables interpolated are p, u,*, p, t^j , and 
e. The min-mod limiter [26] is used to avoid spurious oscillations in the neighbor- 
hood of a discontinuity. Other limiters are available and will be investigated in future 
work to improve monotonicity and convergence behavior. 

The remaining terms to be discretized are the diffusive fluxes and the source 
terms. Consistent with the elliptic nature of the diffusive fluxes, a finite- volume 
representation of a second-order accurate central-difference operator [26, 27] is em- 
ployed. Derivatives required in the diffusive flux evaluation at the cell interface are 
approximated with Gauss’s divergence theorem, integrating around an auxiliary cell 
centered at the interface. Flow variables required at this interface are obtained from 
arithmetic averaging of neighboring cell averages. Derivatives required for the source 
terms are also calculated using Gauss’s divergence theorem by integrating around the 
computational cell. 

To accomodate geometrically complex configurations, we implement a multi-block 
procedure which requires C° grid continuity. 
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4.2 Roe Flux-Difference Splitting 

The interface flux for the finite volume formulation is calculated in each of the three 
coordinate directions as the solution of a locally one-dimensional Riemann problem 
normal to the cell interface (the so-called operator splitting approach) 


1 ^l F =° 

(4.4) 

Q = Ql n < 0 
Q = Qr n > 0 

(4.4a) 


where n is the coordinate normal to the cell interface (£, 77 or £) and 

puU + h x p 

pvU + hyp 

pwU + h z p 

pUH 

pJJ ^~xx 

pUjyy 

pUjzz 

pU T x y 

pUjxz 

piJ Ty Z 

pUc 

U = h x u + h y v + h z w. (4.4c) 

The solution of the Riemann problem results in a shock wave, a contact discontinuity, 
and a rarefaction evolving in time at the interface (see Fig. 1). The interface flux 
can then be determined as the flux at the left or right state incremented by the flux 
differences crossed from that state to the interface. The exact solution of the Riemann 
problem requires an iterative procedure and is quite expensive. A cheaper alternative 
developed by Roe is to construct the solution to the approximate, linearized problem: 

|Q + 4<3 = 0 (4-5) 

where the Jacobian matrix, A = dF jdQ, is evaluated at an average state such that 
it satisfies the jump condition between the flux states at the right and left: 

F r -F l =J(Q r -Ql ) (4.6) 

(The complete form of the Jacobian matrix and the symmetrization matrices is given 
in the Appendix). The Jacobian matrix, A, has twelve eigenvalues, A,-, 

Ai = A2 — A3 = . . . Aj 2 = f/|Vn| 
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A 4 = (u + a) |Vn| (4.7) 

X 5 = (U - a) |Vn|. 

The Roe averaged variables that satisfy Eq. (4.6) are 

p — \J ~PlPr (4.8a) 

fr _ )/PlHl + s/p^Ur ( 4 . 86 ) 

VpZ+VpZ 

= Vplul ± y/fWR (4.8c) 

yfpL + y/PR 

$ = yfnvL_+ ( 4 . 8 rf) 

yfpi + y/^fl 

_ + \fpR w R (4.8c) 

V^+V^ 

q 2 = u 2 + n 2 + w 2 (4-8 /) 

a 2 = (7—1 )(H - q 2 j2 - k) (4-8 g) 

— _ y/PLTjjl, + \fpn T ijH (4 gM 

y/ZZ+y/PZ 

k = {tZ + fQ + f^)/2 (4-8*) 

_ y//>L<T + y/pR^R (4 gj) 

y/pl + y/^H 


The interface flux can finally be written as the average of the interface flux calcu- 
lated from the left state crossing negative running waves and the right state crossing 
positive running waves, as: 

Ti+if 2 = \\F r + F l - |t| ( Qr - Ql )] (4.9) 

This may be written in a more computationally efficient manner as 

F t+ 1/2 = l -[F R + F l - £ |AF|] (4.10) 
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where the flux differences can be derived as: 
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4.3 Time Integration 


The conservation and transport equations given by Eq. (4.1) must be integrated 
in time to provide a steady-state solution. The Reynolds stress and two-equation 
turbulence models are known to exhibit stiff behavior, especially in the near-wall 
region. Therefore, an implicit time integration is used to integrate the equations. 
The Euler implicit time integration scheme can be written in delta form as: 


JAt 


+ S ( (A - Ay) + S V (B - B v ) + S c (C -Cy)-D 


A Q =-R n 


(4.11) 


where A Q = <2" +1 - Q n and A = 8F/8Q, B = dG/dQ, C = dH/dQ , A, = dFJdQ, 
B v = dGyjdQi C v = dH v /dQ, D = dS/dQ are the Jacobian matrices and R n is the 
steady-state residual at time n. Eq. (4.11) results in a large banded system of block 
matrices which must be inverted at each time step to update the solution. 

The storage requirements to solve Eq. (4.11) exactly for a large three-dimensional 
problem are prohibitive. Therefore, an approximate solution at each time step is 
obtained by solving a factored form of the equation. A spatially split, approximate 
factorization scheme can be written as the following four sweeps (the fourth sweep is 
a point implicit treatment of the source terms) through the flowfield: 


[. I - JAtD] A Q' 

= -JAtR n 

[I+JAt6t(A-A v )] A Q” 

= A Q' 

[I + JAtS n {B - B V ))AQ"' 

= A Q" 

[I + JAt6((C — C V )]AQ 

= A Q'" 

Q n+1 

= Q n + AQ 


(4.12) 


The first factor requires only a block inversion at each point in the flowfield. The 
three spatial factors require the inversion of a block pentadiagonal system (block 
12 x 12 matrices for Reynolds stress models or block 7x7 matrices for two-equation 
models) in each of the three coordinate directions (£, 7 /, () as the computational stencil 
for the inviscid fluxes spans five cell-centers. A computational savings can be gained 
by treating the inviscid fluxes as first-order accurate in the implicit sweeps. This 
reduces the system to block tridiagonal matrices. There is no reduction in accuracy at 
steady-state due to the implementation in delta form. A second benefit of first-order 
implicit differencing of the inviscid fluxes is gained from the unconditional diagonal 
dominance of the first-order system. 

The source terms are currently treated in a point implicit manner in the above 
algorithm. Recent work on approximate factorization schemes with source terms [28] 
recommends including the source terms in the spatial factors rather than as a separate 
inversion to reduce factorization error. The current scheme has been chosen to avoid 
storing the source Jacobian and to facilitate the implementation (at a future date) of 
a diagonalized version of the scheme [29]. 

The source term is a function of both the conserved variables and derivatives 
of the conserved variables, S — S(Q,dQ /<?£,). For the current calculations, D is 
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simplified to include only the implicit contributions from the conserved variables in 
order to preserve the point implicit treatment. The inclusion of the dS/d(dQ /dxi) 
terms is possible in the spatially split factors, but is not included at present. The net 
result of neglecting the contribution due to the derivatives of the conserved variables 
is to essentially decouple the turbulence transport equations from the mean flow 
equations in the source term treatment. This is exactly what most researchers have 
implemented when they solve the mean flow equations coupled and then update the 
turbulence tranpsort equations separately. However, there is still a coupling of the 
normal stress components with the mean flow equations as can be seen from the 
Jacobian matrices. 
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5 Results 


As mentioned previously, the intention in this report is not to evaluate turbulent 
closure models but to develop a numerical methodology which is capable of using 
either two-equation or Reynolds stress closure models in the solution of compressible 
flow equations. The task of evaluating closure models in an unambiguous fashion 
requires the use of a common numerical procedure so that any bias from the numerics 
can be consistently minimized. As was discussed in the previous section, there are a 
variety of ways of handling both the coupling of the equations and the differencing 
of the various terms. With the development of the technique outlined in this report, 
it will now be possible to evaluate a variety of closure models. This evaluation is 
presently just beginning and will be the subject of future reports. 

In this section, the results from some initial computations using both two-equation 
and Reynolds stress models will be presented. These are not intended as a compari- 
son of model performance, rather these are intended to illustrate the numerical code’s 
capability. In fact, a fair comparison of closure models in wall-bounded compress- 
ible flows may not be possible at this time. As has been alluded to earlier, the high 
Reynolds number form of both two-equation and Reynolds stress models has been 
fairly well established, even in the compressible regime for moderately high Mach 
numbers. However, the development of appropriate near-wall models for these clo- 
sures in order to integrate directly to the wall has not been completed, particularly in 
the compressible high Mach number regime. This is particularly crucial when one real- 
izes that the stiffness problems commonly associated with the Reynolds stress models 
are due to an improper balance of terms in the modeled equations for the near-wall 
region. Another factor complicating any sensible comparison of two-equation and 
Reynolds stress models is the lack of experimental data. Even in cases where there 
are data, the data are by no means complete, that is, measurements of turbulence 
profiles are rare and without them any comparisons of model performance is indi- 
rectly obtained through comparisons with mean flow variables. The approach taken 
here is to use a variable density extension of an incompressible near-wall model and 
to analyze its ability to calculate the turbulent stress components. The k — e version 
of the Speziale, Abid and Anderson model [6] will be used for the two-equation tests 
and the Shima model [20] will be used for the Reynolds stress test. 

One test case that will be presented is the two-dimensional flow over a 10° com- 
pression ramp. The ramp flow is Mach 3 with a mean flow Reynolds number of 
10 x 10 6 and adiabatic wall conditions. The free-stream temperature, T [», is set at 
300° K. The flow field will be calculated using both the two-equation and Reynolds 
stress models. 

A second test case is the three-dimensional flow over a 5° cone at 2° angle of 
attack. The cone flow is Mach 3.5 with a mean flow Reynolds number of 2.21 x 10 7 
and adiabatic wall conditions. The free-stream temperature is 311°A'. The cone 
flow field will be run with the two-equation model only just to validate the three- 
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dimensional capability of the code. 

The boundary conditions for the calculations are applied explicitly. The inflow 
data are completely set to free-stream values for supersonic calculations. The free- 
stream turbulent kinetic energy level is specified as 1%, and the free-stream turbulent 
eddy viscosity is assumed equal to the free-stream molecular viscosity. From the 
free-stream value of turbulent kinetic energy and turbulent eddy viscosity, the free- 
stream value of solenoidal dissipation rate is obtained from Eq. (2.6a). An isotropic 
distribution is assumed for the normal components of the Reynolds stress tensor 
with the shear stress components set to zero for consistency. Outflow boundary 
conditions are zeroth-order extrapolations of the variables. The farfield boundary 
conditions for mean variables are held to free-stream values since the grid domain 
extends outside the bow shock of the flat-plate leading-edge (for the ramp problem) 
or the cone apex. Turbulence variables are obtained from zeroth-order extrapolation 
from the interior. This Neumann condition is used to minimize any farfield bias to 
the interior turbulence levels. The wall boundary conditions are specified as zero 
velocity (u, = 0), zero Reynolds stresses (r,j = 0), and adiabatic or specified wall 
temperature. At the wall, the pressure is extrapolated from the interior solution with 
a zeroth-order extrapolation, and the solenoidal dissipation rate is equated to the 
second derivative of the turbulent kinetic energy (a consequence of evaluating the 
turbulent kinetic energy equation at the wall). 

The results from the ramp calculations will be shown first. The domain was 
discretized using a 101 by 51 mesh, with high grid clustering in the near-wall region, 
for both the two-equation and Reynolds stress calculations. It does not appear to 
be worthwhile to do grid refinement studies in the context of the present report 
since the closure models have not been finalized nor has a suitable (sufficiently well 
documented) set of test cases been chosen here. These are the topic of future research 
and are presently being initiated. Nevertheless, with the same grid structure, it will 
be possible to get a qualitative comparison of both the two-equation and Reynolds 
stress closures. The mean pressure and velocity, turbulent Reynolds stresses (kinetic 
energy) and turbulent dissipation rate variables that will be shown are all normalized 
by the free-stream mean pressure, free-stream mean density, free-stream mean velocity 
and flat-plate/ramp length. 

Pressure contours for the 10° ramp using the two-equation model are shown in 
Fig. 2. Note that the domain shown is inclusive of the plate leading-edge (lower left 
corner of figure). The shock is captured very cleanly with no apparent oscillations, 
validating the shock capturing properties of the upwinding procedure for the coupled 
set of equations. The apparent thickening of the shock at the outer extent of the 
domain is due to the coarse grid in this region resulting from high grid stretching to 
provide adequate resolution in the near-wall region. The shock appears to affect the 
boundary layer near the start of the ramp and then has minimal direct affect farther 
downstream; however, as subsequent plots will show, the start of the ramp will have 
an affect on the flow field farther downstream along the ramp. 
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The mean velocity profiles are shown in Fig. 3. These profiles are the Favre- 
averaged Cartesian mean velocities — (£f,tq 0), normalized by the free-stream ve- 
locity Uoo, with the domain set between x = 0.0 and x = 1.0. The start of the ramp 
is located at x = 0.5. The four stations show one set of profiles upstream of the ramp 
(a; = 0.396), a second at the start of the ramp (x = 0.499), a third slightly down- 
stream of the ramp (x — 0.598), and a fourth further downstream along the ramp 
( x = 0.808). The figure shows the dominance of the x-component of velocity over the 
entire range, although after the ramp, the y-component does increase significantly 
over its flat-plate value. At the first station, the u-velocity is essentially coincident 
with the vertical axis. The small backflow at the start of the ramp can be seen in the 
x — 0.499 plot. 

Streamwise variation of the turbulent kinetic energy is shown in Fig. 4. The 
downstream development of the turbulent kinetic energy is clearly altered by the 
presence of the ramp. The ramp causes an increase in the turbulent energy away 
from the wall at the expense of the peak energy level near the wall (x = 0.499). 
Further downstream, the profile recovers and tends toward its flat-plate distribution. 

Figure 5 shows the corresponding turbulent dissipation rate evolution. Once again, 
there is seen a significant effect on the turbulence as the flow is subjected to the 
compression ramp. At x = 0.499, the peak in the dissipation rate is shifted away 
from the wall and then by x « 0.808 it has recovered just as in the turbulent kinetic 
energy case. Unfortunately, the quantitative aspects of these results are questionable, 
since the dissipation rate should have its maximum value at the wall (e.g. [10]). 

Clearly, the model used is unable to capture this trend and needs to be improved. 
Such improvements are the focus of research at this time, and as these improvements 
in the models are developed they will be incorporated into the numerical code. 

As a further validation of the code capabilities, the previous ramp flow was also 
computed using the Reynolds stress turbulence model of Shima [20]. The pressure 
contours for the flow are shown in Fig. 6. The results are qualitatively consistent 
with the results from the two-equation model. As was shown in the previous results, 
the shock itself does not affect the turbulence quantities within the boundary layer 
flow, but the effects of the ramp start will have an effect on the turbulence quantities 
both at the start as well as farther downstream along the ramp. Again the shock is 
captured cleanly with minimal oscillations, verifying the shock capturing capability 
of the scheme including Reynolds stress models. The shock thickening at the outer 
region of the domain is again due to grid coarseness in this region. 

Mean velocity profiles are shown in Fig. 7 at the same streamwise stations shown 
for the two-equation model. The start of the ramp has the effect of retarding the flow 
near the wall, with the u velocity profile at x = 0.499 showing that the flow is near 
incipient separation. The two-equation results did show separation at this point, but 
it was rather weak. Clearly, the initiation of the ramp causes a strong deceleration of 
the flow and this affects the turbulence at the ramp start as well as downstream. The 
velocity profiles computed using the Reynolds stress model are fuller than those from 
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the two-equation model suggesting that the predicted wall skin friction is higher. As 
in Fig. 3, the u-component velocity only begins to have an impact at the ramp start, 
and in this case the effect on the flat-plate and near the ramp start are minimal. 

The turbulent kinetic energy shown in Fig. 8 is obtained from the trace of the 
normal stress components of the Reynolds stress tensor. The qualitative features of 
the profiles are similar to those obtained using the two-equation model. The effect of 
the ramp once again causes the turbulent kinetic energy to increase farther from the 
wall and then to subsequently relax back farther downstream. A comparison with the 
two-equation results in Fig. 4 clearly shows that the Reynolds stress closure produces 
a significantly smaller response to the ramp start than the two-equation closure. The 
difference in magnitude between the two-equation and Reynolds stress results is easily 
explained. In the two-equation case, the Reynolds stress components are extracted 
from the Boussinesq approximation which causes a ‘rapid response’ of the stresses 
to the abrupt change of geometry at the ramp start; whereas, the Reynolds stress 
model has inherent relaxation and memory characteristics that cause a more tempered 
response to such geometric changes. 

Figure 9 shows the turbulent dissipation rate profiles at various streamwise sta- 
tions. These results suffer from the same deficiency that the dissipation rate from the 
two-equation model had in that the dissipation rate does not peak at the wall. Other- 
wise, the trends are the same as in Fig. 5 where the ramp caused the dissipation rate 
to increase initially and then relax back to its flat-plate level farther downstream. In 
this case, the turbulent dissipation rate levels are less than the corresponding levels 
for the two-eqaution closure. This is, of course, consistent with trends shown for the 
turbulent kinetic energy. 

An advantage of the Reynolds stress closure model is the ability to give the be- 
havior of the component stresses. The normal stress components are shown in Fig. 10 
where the t xx , r yy and t zz Reynolds stresses are shown. As the figure shows, the t xx 
and r zz components dominate in their contribution to the turbulent kinetic energy. 
The start of the ramp initiates a large increase in the level of component energy from 
the flat-plate levels, which as the kinetic energy plot shows, persists dowmstream to 
the x = 0.598 station. All three components show an increase in energy level over 
a wide region of the boundary layer at this station. This ‘overshoot has diminished 
somewhat by the station x = 0.808. 

The only component of the Reynolds shear stress to survive in this two-dimensional 
flow is the T xy component. Figure 11 shows its variation with downstream distance. 
The shear stress profile is also affected by the ramp start, although not to the same 
extent as the other turbulence profiles. In addition, downstream of the ramp start, 
along the ramp, the Reynolds shear stress goes positive in the near- wall region. While 
some regions of positive shear stress may be expected near the point of incipient sep- 
aration, the wide streamwise extent of the positive shear stress values is questionable. 
This widespread effect is probably due to an error in the near- wall model and/or 
the grid resolution in this region. Shima’s Reynolds stress model was developed for 
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incompressible flows and may prove to be inadequate for compressible, high Mach 
number flow fields with a simple variable density extension. Downstream of the ramp 
start, the shear stress levels relax back to their flat-plate levels. 

The next set of results to be shown are for the three-dimensional turbulent flow 
over an axisymmetric cone at angle of attack. The domain was discretized into a 49 
by 65 by 25 mesh (streamwise by normal by circumferential), with high grid clustering 
once again near the solid surface. Three-dimensional calculations are computationally 
intensive and for the purposes of the present report, only a two-equation turbulent 
closure model was applied to this flow. The variables shown will be normalized as 
before in the ramp case. 

Figure 12 shows the cone at angle of attack to the mean flow with pressure contours 
superimposed. The bow shock is again captured with no apparent oscillations. The 
grid has extreme resolution in the near wall region (approximately 20 grid points 
are within the region y + < 20), resulting in an extremely coarse grid at the outer 
boundary and thick bow shock. Inviscidly, this cone flow field would exhibit conical 
scaling, where the solution would be constant along rays eminating from the cone 
apex. Due to viscous effects, the flow field exhibits only quasi-conical behavior away 
from the cone apex, allowing us to look at a single cross-section of the calculation to 
gain qualitative, but not quantitative, information about the entire flow field. For this 
purpose, the station x = 0.83 (based cone length of unity), which is well downstream 
of the cone apex, is chosen to show the azimuthal variation of the mean and turbulent 
quantities. 

The Cartesian mean velocity components ft,- = ( u , v, w) (scaled by the free-stream 
velocity) are shown in Fig. 13 at four azimuthal locations: (j> = 86.25° (top, leeward 
side), 4> = 26.25°, <j> = -26.25°, and <f> = -86.25° (bottom, windward side). The ver- 
tical coordinate is the radial distance measured from the x-axis. In all the figures, the 
normal and circumferential velocities have minimal effect on the total mean velocity. 
From the streamwise component profile, it appears that there is a strong suppression 
of the turbulence on the windward side of the cone where the boundary layer is quite 
thin. On the leeward side, a thicker boundary- layer region is present, since the layer 
is not being constrained by a strong impinging mean flow. 

Figure 14 shows the corresponding azimuthal variation of the turbulent kinetic 
energy. The figure shows that the boundary layer is quite thin on both the leeward 
and windward sides of the cone, although the leeward side boundary layer thickness 
is about twice the size of the windward side. The peak turbulent intensity, which 
occurs very near the wall, is relatively constant over the azimuthal range of the plot. 

The solenoidal turbulent dissipation rate is shown in Fig. 15. As expected, the 
trends are consistent with the kinetic energy profiles, with the peak dissipation rate 
increasing as the windward side of the cone is approached. Note that significant 
dissipation levels are constrained to regions very near the wall. 
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6 Concluding Remarks 

The purpose of this report is to establish both the differential and numerical 
framework for the development of a numerical algorithm for the solution of compress- 
ible turbulent flow problems with both two-equation and Reynolds stress turbulence 
closure models. As the report has shown, the task is straightforward, but rather cum- 
bersome. There are a significant number of unknown correlations which need to be 
modeled and which have not received sufficient attention by the scientific community 
to be adequately validated. This validation process, however, can only be accom- 
plished by the development of the type of numerical code outlined here. With the 
present tool, a variety of complex flows can be examined using both the two-equation 
and Reynolds stress closure models with the dual goal of developing better closure 
models and comparing the predictive capabilities of both two-equation and Reynolds 
stress formulations. 

The test problems that were computed here were chosen to show the capabilities 
of the code rather than to compare the models used. If a formal comparison were 
intended, improved closure models would need to be used. Nevertheless, even with 
this type of qualitative study, results were obtained for both the ramp and cone 
problems which showed some interesting chararcteristic features of the mean and 
turbulent flow field. Unfortunately, a more thorough study is rather computationally 
costly and should only be performed after the methodology is validated. Work is 
beginning on this next phase of research where improved versions of the two-equation 
and Reynolds stress models will be implemented into the code which will be used for 
predicting selected flow problems. 
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Appendix 

The Jacobian matrices, (.4, 5,(7), may be written for n = (£,??,£) as: 
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The Jacobian matrix may be written as = TAT h The symmetry matrix, T, 
iven as: 


is given as: 
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where 


V = l X U + lyV + l z W 
TT = rh x u + rhyV + 
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The inverse of the symmetry matrix, T x , is given as: 
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where 

T = ( 7 -l)/a 2 

The vectors l and m are non-unique, mutually orthogonal to n, and tangent vec- 
tors in the plane of the cell interface. The flux differences given in Section 4 use 
metric identities to provide a form which is independent of these tangency vectors for 
computational efficiency. 
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Figure 1: Diagram of shock and rarefaction wave solution to Riemann problem. 




Figure 2: Pressure contours for Mach 3, 10° compression ramp using a two-equation 
turbulence model. 




V.LTUV U. l/UU 

0.0 0.4 0.8 0.0 0.4 0.8 0.0 0.4 0.8 0.0 0.4 0.8 

Mean Velocities 


Figure 3: Streamwise variation of mean velocity profiles for 10° compression ramp 
using a two-equation turbulence model: u/U 0 c , ; v/Uoo , — — . 
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Turbulent Kinetic Energy (*10 2 ) 

Figure 4: Streamwise variation of turbulent kinetic energy profiles for 10° compression 
ramp using a two-equation turbulence model. 
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Turbulent Dissipation Rate 

Figure 5: Streamwise variation of turbulent dissipation rate profiles for 10° compres- 
sion ramp using a two-equation turbulence model. 
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Turbulent Kinetic Energy (xIO 2 ) 

Figure 8: Streamwise variation of turbulent kinetic energy profiles for 10° compression 
ramp using a Reynolds stress turbulence model. 
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Figure 9: Streamwise variation of turbulent dissipation rate profiles for 10° compres- 
sion ramp using a Reynolds stress turbulence model. 
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Figure 10: Normal Reynolds stress component profiles for 10° compression ramp. 
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Figure 11: Reynolds shear stress profiles for 10° compression ramp. 
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Figure 12: Pressure contours for Mach 3.5, 5° cone at 2° angle of attack using a 
two-equation turbulence model. 
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Figure 13: Azimuthal variation of mean velocities for 5° cone at 2° angle of attack 
using a two-equation turbulence model: ujUoo , ; v/Uoo, — — ; w/Uoo, - - -. 
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Turbulent Kinetic Energy (xIO 2 ) 

Figure 14: Azimuthal variation of turbulent kinetic energy for 5° cone at 2° angle of 
attack using a two-equation turbulence model. 
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Figure 15: Azimuthal variation of turbulent dissipation rate for 5° cone at 2° angle 
of attack using a two-equation turbulence model. 
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